#include "Modeling\CUDADataStruct.h"
//#include ".\Modeling\RenderBaseStruct.h"
#include ".\CUDA\utility\inc\cutil.h"
#include ".\CUDA\utility\inc\cudpp\cudpp.h"
#include <stdio.h>

typedef unsigned int uint;

#define UNIFORM_KNOT 1
#define SORT_SHARE   1
#define NORMALFUNCTION	0

#define NTHREADS		128
#define PIXELTHREADS	64
#define DEBUG_BUFFERS 0

#define ERROR_FUNC	1.0e-6
#define ERROR_DELTA	1.0e-3

extern void dumpBuffer(uint *d_buffer, int nelements);
extern void dumpBufferfloat(float *d_buffer, int nelements);

const uint pixelBits = 10;
uint patchBits, pixelMask, patchMask;


uint3 deg, degShift, seg, segShift;
uint degSize, segSize;
float3 boxSize, box0;

uint nPatch;

#ifndef UNIFORM_KNOT
__constant__ float* dBoxX = 0;
__constant__ float* dBoxY = 0;
__constant__ float* dBoxZ = 0;
#endif

__constant__ int* dPatchIndex = 0;
__constant__ int* dPatchIndexInverse = 0;
__constant__ float* dFuncData = 0;

__device__   int* dPolarIndex = 0;
__device__	 float* dPolarData = 0;
texture<float, 1, cudaReadModeElementType> polarTex;


//the variable below are defined in the way described in "Progress 06.22"
uint maxSeg;
uint maxSilPoint;
uint actualSilPoint;

uint allPointSize;
uint maxPoint;
uint actualPoint;


__device__ size_t* dMaxSilPoint  = 0;
__device__ size_t* dMaxPoint = 0;

__device__ uint* dSilOccupied = 0;
__device__ uint* dPointOccupied = 0;
__device__ uint* dSilOccupiedValid = 0;
__device__ uint* dPointOccupiedValid = 0;
__device__ uint* dOccupiedCompact = 0;
__device__ uint* dOccupiedSort = 0;
__device__ uint* dOccupiedSortPad = 0;

uint* pointOccupiedSort = 0;
uint* pointOccupiedPad = 0;

CUDPPHandle plan;

__constant__ MULTI_SIL_RENDER* dMultiSilRender;
int renderWidth, renderHeight;
//float4 viewProj0, viewProj1, viewProj2, viewProj3;
float4 invViewProj0, invViewProj1, invViewProj2, invViewProj3;
//float3 eye;
//float errorXTol;
//LIGHT light[3];
//MATERIAL material;
//int nLight;

__device__ void normalize(float* a)
{
	float sum;
	sum=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
	if(sum > 1.0e-10f)
    {
		a[0]/=sum; a[1]/=sum; a[2]/=sum;
	}
}

__device__ void addambient(MATERIAL *mat, float rgb[])
{
	rgb[0]+=mat->ka*mat->d[0]*mat->a[0];
	rgb[1]+=mat->ka*mat->d[1]*mat->a[1];
	rgb[2]+=mat->ka*mat->d[2]*mat->a[2];
}

__device__ void addPhong(MATERIAL *mat, LIGHT* light, float norm[], float ray[], float rgb[])
{
	float H[3] ={(-light->dir[0]-ray[0])/2.0f, (-light->dir[1]-ray[1])/2.0f, (-light->dir[2]-ray[2])/2.0f};
	normalize(H);
	float phong = norm[0]*H[0]+norm[1]*H[1]+norm[2]*H[2];
	phong *= phong;
    phong *= phong;
    phong *= mat->ks;
	rgb[0]+=phong*light->rgb[0];
	rgb[1]+=phong*light->rgb[1];
	rgb[2]+=phong*light->rgb[2];
}

__device__ void addDiffuse(MATERIAL *mat, LIGHT* light, float diff, float norm[], float rgb[])
{
	rgb[0]+=mat->kd*diff*light->rgb[0];
	rgb[1]+=mat->kd*diff*light->rgb[1];
	rgb[2]+=mat->kd*diff*light->rgb[2];
}

__device__ void GaussEliminationGPU(float A[][3], float b[])
{
	int i,j,k,ik;
	float mik, temp;
	for(k=0; k<3; k++)
	{
		mik=-1.0f;
		for(i=k; i<3; i++)
		{
			if(fabs(A[i][k])>mik)
			{
				mik=fabs(A[i][k]);
				ik=i;
			}
		}
		if (mik<1.0e-10f)
		{
			b[0]=b[1]=b[2]=0.0f;
			return;
		}
		for(j=k; j<3; j++)
		{
			//swap(A[ik][j], A[k][j]);
			temp = A[ik][j];
			A[ik][j]=A[k][j];
			A[k][j]=temp;
		}
		//swap(b[ik], b[k]);
		temp = b[ik];
		b[ik] = b[k];
		b[k] = temp;
		b[k] /= A[k][k];
		for(i=3-1; i>=k; i--)
			A[k][i]/=A[k][k];
		for(i=k+1; i<3; i++)
		{
			b[i] -= A[i][k]*b[k];
			for (j=3-1; j>k; j--)
				A[i][j] -= A[i][k]*A[k][j];
		}
	}
	for(i=3-1; i>=0; i--)
	{
		for (j=i+1; j<3; j++)
			b[i] -= A[i][j]*b[j];
	}
}

__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}

__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossomPolarTex(float u, float v, float w, int nx, int ny, int nz, int index, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += tex1Dfetch(polarTex, index++) * nn[k];
	}

	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);
	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void EvalBezierBlossomNormal(float u, float v, float w, int nx, int ny, int nz, uint3 degShift, float* func, float* re)
{
	float nn[4];
	float f1[4][4];
	float f2[4];

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += *func++ * nn[k];
	}
	func += degShift.x;
	GetBernstein(nn, ny, v);
	for (int i=0; i<nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx-1, u);
	re[0] = 0.0f;
	for (int i=0; i<nx; i++)
	{
		re[0] += f2[i]*nn[i];
	}

	//X direction:
	GetBernstein(nn, nz, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<=nz; k++)
				f1[i][j] += *func++ * nn[k];
		}
		func += degShift.y;
	}
	GetBernstein(nn, ny-1, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[1] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[1] += f2[i]*nn[i];
	}

	//Z direction:
	GetBernstein(nn, nz-1, w);
	for (int i=0; i<=nx; i++)
	{
		for (int j=0; j<=ny; j++)
		{
			f1[i][j] = 0.0f;
			for (int k=0; k<nz; k++)
				f1[i][j] += *func++ * nn[k];
			func++;
		}
	}
	GetBernstein(nn, ny, v);
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}
	GetBernstein(nn, nx, u);
	re[2] = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		re[2] += f2[i]*nn[i];
	}
}



extern "C"
void InitMulitSilSurfaceCUDA(float* pFuncData, int* Index, int* indexInverse, int NumPatch, MULTI_SIL_SURFACE* multiSilSurface, int surfWidth, int surfHeight)
{
	patchBits = 8*sizeof(uint)-2*pixelBits;
	uint maxNumPatch = 1<<patchBits;
	if (NumPatch>maxNumPatch)
	{
		printf("Too many patchs!\n"); 
		return;
	}

	pixelMask = (1<<pixelBits)-1;
	patchMask = (1<<patchBits)-1;


	nPatch = NumPatch;
	deg.x=multiSilSurface->degX; deg.y=multiSilSurface->degY; deg.z=multiSilSurface->degZ;
	seg.x=multiSilSurface->segX; seg.y=multiSilSurface->segY; seg.z=multiSilSurface->segZ;
	degShift.x = (deg.y+1)*(deg.z+1); degShift.y = deg.z+1; degShift.z = 1;
	segShift.x = seg.y*seg.z; segShift.y = seg.z; segShift.z = 1;
	maxSeg = max(seg.x, seg.y); maxSeg = max(maxSeg, seg.z);
	degSize = (deg.x+1)*degShift.x; segSize = seg.x*segShift.x;

	int funcSize = nPatch*degSize*sizeof(float);	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dFuncData, funcSize) );
	CUDA_SAFE_CALL( cudaMemcpy(dFuncData, pFuncData, funcSize, cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dFuncData, pFuncData, size);	

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPatchIndex, segSize*sizeof(int)) );
	CUDA_SAFE_CALL( cudaMemcpy(dPatchIndex, Index, segSize*sizeof(int), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dPatchIndex, Index, segX*segY*segZ*sizeof(uint));

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPatchIndexInverse, nPatch*sizeof(int)) ) ;
	CUDA_SAFE_CALL( cudaMemcpy(dPatchIndexInverse, indexInverse, nPatch*sizeof(int), cudaMemcpyHostToDevice) );

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPolarData, funcSize) );
	//cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
	CUDA_SAFE_CALL( cudaBindTexture(0, polarTex, dPolarData, cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat)/*channelDesc*/) );
	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dPolarIndex, nPatch*sizeof(int)) );
	
	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dMultiSilRender, sizeof(MULTI_SIL_RENDER)) );
	
	renderWidth = surfWidth;
	renderHeight = surfHeight;
 	allPointSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupied, allPointSize) );		
	CUDA_SAFE_CALL( cudaMalloc((void**)&dSilOccupiedValid, allPointSize) );

	CUDA_SAFE_CALL( cudaMalloc((void**)&dPointOccupied, allPointSize) );		
	CUDA_SAFE_CALL( cudaMalloc((void**)&dPointOccupiedValid, allPointSize) );

	CUDA_SAFE_CALL( cudaMalloc((void**)&dOccupiedCompact, allPointSize) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dOccupiedSort, allPointSize) );

    CUDA_SAFE_CALL( cudaMalloc( (void**) &dMaxSilPoint, sizeof(size_t)));
    CUDA_SAFE_CALL( cudaMalloc( (void**) &dMaxPoint, sizeof(size_t)));

	uint numSortPad = allPointSize+sizeof(uint);
	pointOccupiedSort = (uint*)malloc(numSortPad);
	pointOccupiedPad = (uint*)malloc(numSortPad*5/3);
	printf("Attention: We got memory leak here!\n");
	CUDA_SAFE_CALL( cudaMalloc((void**)&dOccupiedSortPad, numSortPad*5/3) );


#if	UNIFORM_KNOT	
	box0.x = multiSilSurface->boxX[0];
	box0.y = multiSilSurface->boxY[0];
	box0.z = multiSilSurface->boxZ[0];
	boxSize.x = multiSilSurface->boxX[1]-multiSilSurface->boxX[0];
	boxSize.y = multiSilSurface->boxY[1]-multiSilSurface->boxY[0];
	boxSize.z = multiSilSurface->boxZ[1]-multiSilSurface->boxZ[0];
#else	
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxX, (seg.x+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxY, (seg.y+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMalloc((void**)&dBoxZ, (seg.z+1)*sizeof(float)) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float), cudaMemcpyHostToDevice) );	
	CUDA_SAFE_CALL( cudaMemcpy(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float), cudaMemcpyHostToDevice) );
	CUDA_SAFE_CALL( cudaMemcpy(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float), cudaMemcpyHostToDevice) );
	//cudaMemcpyToSymbol(dBoxX, multiSilSurface->boxX, (seg.x+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxY, multiSilSurface->boxY, (seg.y+1)*sizeof(float));
	//cudaMemcpyToSymbol(dBoxZ, multiSilSurface->boxZ, (seg.z+1)*sizeof(float));
#endif
}


//__global__ void GetPolarData(float* FuncData, float* PolarData, int* PatchIndex, int* PolarIndex, float3 eye, uint3 deg, uint3 degShift, uint3 segShift, uint funcSize, float3 box0, float3 boxSize)
__global__ void GetPolarData(float* FuncData, float* PolarData, /*int* PatchIndex, */int* PatchIndexInverse, int* PolarIndex, float3 eye, uint3 deg, uint3 degShift, uint3 segShift, uint funcSize, float3 box0, float3 boxSize)
{
	//int patchIndex = PatchIndex[blockIdx.x];
	//if(patchIndex<0)		return;
//int patchIndex = PatchIndexInverse[blockIdx.x];


	extern __shared__ float funcPolar[];
	float* func = funcPolar;
	float* polar = &func[funcSize];
	float* temp = &polar[funcSize];

	//__shared__ float func[27];
	//__shared__ float polar[27];
	//__shared__ float temp[27];

	float A[4][3] = {{1,0,0}, {0.5f, 0.5f, 0.0f}, {0.0f,1.0f, 0.0f}, {0,0,0}};
	if(deg.x==3&&deg.y==3&&deg.z==3)
	{
		//float A[4][3] = {{1,0,0}, {1.0f/3.0f,2.0f/3.0f,0.0f}, {0.0f,2.0f/3.0f,1.0f/3.0f}, {0,0,1}};
		A[1][0] = 1.0f/3.0f;	A[1][1] = 2.0f/3.0f;
		A[2][1] = 2.0f/3.0f;	A[2][2] = 1.0f/3.0f;
		A[3][2]	= 1.0f;
	}

	uint ind = threadIdx.x;	
	uint polarIndex = PatchIndexInverse[blockIdx.x];
	func[ind] = FuncData[blockIdx.x*funcSize+ind];
	polar[ind] = 0.0f;
	temp[ind] = 0.0f;
	__syncthreads();
		
	/**/
	//X direction
	uint patchIndexX = polarIndex/segShift.x;
	float patchX = patchIndexX*boxSize.x + box0.x;
	float eyeU = (eye.x-patchX)/boxSize.x;
	uint indexX = ind/degShift.x;
	if(indexX<deg.x)
		temp[ind] = deg.x*(eyeU*func[ind+degShift.x] + (1-eyeU)*func[ind]);
	__syncthreads();
	if(deg.x==2&&deg.y==2&&deg.z==2)
		polar[ind]  += A[indexX][0]*temp[ind-indexX*degShift.x] + A[indexX][1]*temp[ind-indexX*degShift.x + degShift.x]; //+A[indexX][2]*polarJacobi[...+2*degShift.x]
	else if(deg.x==3&&deg.y==3&&deg.z==3)
		polar[ind]  += A[indexX][0]*temp[ind-indexX*degShift.x] + A[indexX][1]*temp[ind-indexX*degShift.x + degShift.x] + A[indexX][2]*temp[ind-indexX*degShift.x + 2*degShift.x];

	
	//Y direction
	uint patchIndexY = polarIndex%segShift.x/segShift.y;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float eyeV = (eye.y-patchY)/boxSize.y;
	uint indexY = ind%degShift.x/degShift.y;
	if(indexY<deg.y)
		temp[ind] = deg.y*(eyeV*func[ind+degShift.y] + (1-eyeV)*func[ind]);
	__syncthreads();
	if(deg.x==2&&deg.y==2&&deg.z==2)
		polar[ind] += A[indexY][0]*temp[ind-indexY*degShift.y] + A[indexY][1]*temp[ind-indexY*degShift.y + degShift.y]; //+A[indexY][2]*polarJacobi[...+2*degShift.y]
	else if(deg.x==3&&deg.y==3&&deg.z==3)
		polar[ind] += A[indexY][0]*temp[ind-indexY*degShift.y] + A[indexY][1]*temp[ind-indexY*degShift.y + degShift.y] + A[indexY][2]*temp[ind-indexY*degShift.y + 2*degShift.y];

	
	//Z direction
	uint patchIndexZ = polarIndex%segShift.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;
	float eyeW = (eye.z-patchZ)/boxSize.z;
	uint indexZ = ind%degShift.y;
	if(indexZ<deg.z)
		temp[ind] = deg.z*(eyeW*func[ind+degShift.z] + (1-eyeW)*func[ind]);
	__syncthreads();
	if(deg.x==2&&deg.y==2&&deg.z==2)
		polar[ind] += A[indexZ][0]*temp[ind-indexZ*degShift.z] + A[indexZ][1]*temp[ind-indexZ*degShift.z + degShift.z]; //+A[indexZ][2]*polarJacobi[...+2*degShift.z]
	else if(deg.x==3&&deg.y==3&&deg.z==3)
		polar[ind] += A[indexZ][0]*temp[ind-indexZ*degShift.z] + A[indexZ][1]*temp[ind-indexZ*degShift.z + degShift.z] + A[indexZ][2]*temp[ind-indexZ*degShift.z + 2*degShift.z];

	//Check the existance of the polar surface using the sign variation of the coefficients.
	if(ind<funcSize-1 && polar[ind]*polar[ind+1]<0)
	{
		PolarIndex[blockIdx.x] = polarIndex;
	}		

	//Here, d=16 is suit for case of 3*3*3 and below.
	//3*3*4 and above should use d=32;
	//temp[ind] = fabsf(polar[ind]);
	//temp[funcSize+ind] = 0.0f;
	//for(int d=16; d>0; d >>= 1)
	//{
	//	__syncthreads();
	//	if(ind<d)
	//		temp[ind] = max(temp[ind], temp[ind+d]);
	//}
	//__syncthreads();
	//PolarData[blockIdx.x*funcSize+ind] = polar[ind]/temp[0];
	PolarData[blockIdx.x*funcSize+ind] = polar[ind];


}

__global__ void NoSortCheckSilPointUniform(unsigned char* surface, size_t pitch, int width, int height, uint pixelBits, uint* dSilOccupied, uint* dSilOccupiedValid)
{	
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	int y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=width || y>=height)		return;
	uint pixelI = y*width+x;
    float* pixel = (float*)(surface + y*pitch) + 4*x;
	if( pixel[0] == 0.0f || pixel[0] == 1.0f)		return;

	dSilOccupied[pixelI] = ((x<<pixelBits)|y);
	dSilOccupiedValid[pixelI] = 1;
}

__global__ void CheckAllSilPointUniform(unsigned char* surface, size_t pitch, int width, int height, uint pixelBits, uint* dSilOccupied, uint* dSilOccupiedValid,
										float3 box0, float3 boxSize, uint3 seg, uint3 segShift, uint patchBits, int* dPatchIndex, uint* dPointOccupied, uint* dPointOccupiedValid,
										float4 invViewProj0, float4 invViewProj1, float4 invViewProj2, float4 invViewProj3)
{	
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	int y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=width || y>=height)		return;
	uint pixelI = y*width+x;
    float* pixel = (float*)(surface + y*pitch) + 4*x;
	//if( pixel[0] == 0.0f || pixel[0] == 1.0f)		return;
	if( pixel[0]<0.1f)								return;

	if( pixel[0]==1.0f )
	{
		//dSilOccupied[pixelI] = ((x<<pixelBits)|y);
		//dSilOccupiedValid[pixelI] = 1;

		float posx = 2.0f*(x+0.5f)/width-1.0f;
		float posy = 1.0f-2.0f*(y+0.5f)/height;
		float posz = pixel[2];
		float worldW =  posx*invViewProj0.w + posy*invViewProj1.w + posz*invViewProj2.w + invViewProj3.w;
		float worldX = (posx*invViewProj0.x + posy*invViewProj1.x + posz*invViewProj2.x + invViewProj3.x) / worldW;
		float worldY = (posx*invViewProj0.y + posy*invViewProj1.y + posz*invViewProj2.y + invViewProj3.y) / worldW;
		float worldZ = (posx*invViewProj0.z + posy*invViewProj1.z + posz*invViewProj2.z + invViewProj3.z) / worldW;

		worldX = (worldX-box0.x)/boxSize.x;
		worldY = (worldY-box0.y)/boxSize.y;
		worldZ = (worldZ-box0.z)/boxSize.z;

		if(worldX>=0 && worldX<seg.x && worldY>=0 && worldY<seg.y && worldZ>=0 && worldZ<seg.z)
		{
			//int patchI = dPatchIndex[(int)floorf(worldX)*segShift.x+(int)floorf(worldY)*segShift.y+(int)floorf(worldZ)*segShift.z];
			int patchI = (int)floorf(worldX)*segShift.x+(int)floorf(worldY)*segShift.y+(int)floorf(worldZ)*segShift.z;
			//if(patchI>0)
			{
				dPointOccupied[pixelI] = (((x<<pixelBits)|y)<<patchBits)|patchI;
				dPointOccupiedValid[pixelI] = 1;
			}
		}
		return;
	}

	dSilOccupied[pixelI] = ((x<<pixelBits)|y);
	dSilOccupiedValid[pixelI] = 1;
}

__global__ void NoSortUpdateSilPointKernel(
											uint* dPointOccupied, uint* dPointOccupiedValid, uint patchBits,
											unsigned char*surfaceIn, unsigned char*surfaceOut, size_t pitch,
											MULTI_SIL_RENDER* dMultiSilRender, uint pixelBits,
											uint* OccupiedCompact, float* FuncData, int* PatchIndex, int* PolarIndex,
											uint3 deg, uint degSize, uint3 seg, uint3 segShift,
											float3 boxSize, float3 box0//,float* elimation
										   )
{
	uint pixelX = OccupiedCompact[threadIdx.x + blockIdx.x*blockDim.x];
	uint pixelY = pixelX<<(sizeof(uint)*8-pixelBits)>>(sizeof(uint)*8-pixelBits);
	pixelX = pixelX>>pixelBits;
	float* pixelInfo = ((float*)(surfaceIn+pixelY*pitch) + pixelX*4 +2);

	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *(pixelInfo);

	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float worldW;
	int patchIndexX, patchIndexY, patchIndexZ;
	float patchX, patchY, patchZ;



	int patchIndex;
	float errF = 1000.0f;
	float deltaPoint = 1000.0f;
	float delta = 0.01f;
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	int step[3] = {0, 0, 0};
	
	int pixelRefine, pixelIter;
	int dir = 0; 
	int dirMin = 0;
	float pad;
	float* pixel;
	
	float plane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	int polarIndex;
	bool isOK = true;

	worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	patchX = patchIndexX*boxSize.x + box0.x;
	patchY = patchIndexY*boxSize.y + box0.y;
	patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;


/**/
	//if(pixelInfo.x<0.3f)
	{
		plane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
		plane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
		plane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
		plane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	}
	//else
	//{
	//	plane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	//	plane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	//	plane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	//	plane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	//}
	plane[0] *= boxSize.x;
	plane[1] *= boxSize.y;
	plane[2] *= boxSize.z;

	//OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = patchIndexX;
	//return;

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];
		
		step[0] = floorf(xyz[0]);
		step[1] = floorf(xyz[1]);
		step[2] = floorf(xyz[2]);

		patchIndexX += step[0];
		patchIndexY += step[1];
		patchIndexZ += step[2];

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		polarIndex = PolarIndex[patchIndex];
		if(polarIndex<0)		break;
		{
			//dir = ((xyz[0]-step[0])>(xyz[1]-step[1]))?0:1;
			//dirMin = 1-dir;
			//dir = ((xyz[dir]-step[dir])>(xyz[2]-step[2]))?dir:2;
			//dirMin = ((xyz[dirMin]-step[dirMin])<(xyz[2]-step[2]))?dirMin:2;			
			//
			//if(xyz[dirMin]-step[dirMin]<0.05f)
			//{
			//	patchIndexX -= step[0];
			//	patchIndexY -= step[1];
			//	patchIndexZ -= step[2];
			//	fValue[0] = fValue[1] = fValue[2] = 0.0f;
			//	fValue[dirMin] = 0.15f;
			//	continue;
			//}

			//if(xyz[dir]-step[dir]>0.95f)
			//{
			//	patchIndexX -= step[0];
			//	patchIndexY -= step[1];
			//	patchIndexZ -= step[2];
			//	fValue[0] = fValue[1] = fValue[2] = 0.0f;
			//	fValue[dir] = -0.15f;
			//	continue;
			//}	
			
			//break;
		}
		patchIndex *= degSize;
		
		patchX += step[0]*boxSize.x;
		patchY += step[1]*boxSize.y;
		patchZ += step[2]*boxSize.z;

		plane[3] += step[0]*plane[0] + step[1]*plane[1] + step[2]*plane[2];

		xyz[0] -= step[0];
		xyz[1] -= step[1];
		xyz[2] -= step[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &fValue[1]);
		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];

		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		EvalBezierBlossomPolarTex(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][1]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, patchIndex, &jacobi[1][2]);

		jacobi[1][0] = (jacobi[1][0]-fValue[1])/delta;
		jacobi[1][1] = (jacobi[1][1]-fValue[1])/delta;
		jacobi[1][2] = (jacobi[1][2]-fValue[1])/delta;

		jacobi[2][0] = plane[0];
		jacobi[2][1] = plane[1];
		jacobi[2][2] = plane[2];


		//elimation[eliIndex  ] = jacobi[0][0];
		//elimation[eliIndex+1] = jacobi[0][1];
		//elimation[eliIndex+2] = jacobi[0][2];
		//elimation[eliIndex+3] = jacobi[1][0];
		//elimation[eliIndex+4] = jacobi[1][1];
		//elimation[eliIndex+5] = jacobi[1][2];
		//elimation[eliIndex+6] = plane[0];
		//elimation[eliIndex+7] = plane[1];
		//elimation[eliIndex+8] = plane[2];
		//elimation[eliIndex+9] = fValue[0];
		//elimation[eliIndex+10] = fValue[1];
		//elimation[eliIndex+11] = fValue[2];
		//patchIndex+=15;
		//for(int i=0; i<15; i++)
		//	elimation[eliIndex+i] = tex1Dfetch(polarTex, patchIndex++);

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];

		//if(its==4)
		//{
		//OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = int(fValue[1]);
		//return;
		//}

		if (errF <= ERROR_FUNC)					break;

		GaussEliminationGPU(jacobi, fValue);

		//elimation[eliIndex+12] = fValue[0];
		//elimation[eliIndex+13] = fValue[1];
		//elimation[eliIndex+14] = fValue[2];


		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= ERROR_DELTA)			break;	
	}

	//if(pixelY==62)
	//{
	//	OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = patchIndexX;
	//	return;
	//}
	//else
	//{
	//	OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 0;
	//	return;
	//}

	if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 0;
		isOK = false;
		//return;
	}

	if(patchIndex<0)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1;
		isOK = false;
		//return;
	}
		
	if(polarIndex<0)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 2;
		isOK = false;
		//return;
	}

	if(errF>ERROR_FUNC)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 3;
		isOK = false;
		//return;
	}

	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	//if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*1000 && deltaPoint<dMultiSilRender->errorXTol*10) )
	if(isOK)
	{
		xyz[0] = xyz[0]*boxSize.x+patchX;
		xyz[1] = xyz[1]*boxSize.y+patchY;
		xyz[2] = xyz[2]*boxSize.z+patchZ;

		fValue[0] = xyz[0]*dMultiSilRender->viewProj[0][0] + xyz[1]*dMultiSilRender->viewProj[1][0] + xyz[2]*dMultiSilRender->viewProj[2][0] + dMultiSilRender->viewProj[3][0]; 
		//fValue[1] = xyz[0]*dMultiSilRender->viewProj[0][1] + xyz[1]*dMultiSilRender->viewProj[1][1] + xyz[2]*dMultiSilRender->viewProj[2][1] + dMultiSilRender->viewProj[3][1]; 
		fValue[2] = xyz[0]*dMultiSilRender->viewProj[0][3] + xyz[1]*dMultiSilRender->viewProj[1][3] + xyz[2]*dMultiSilRender->viewProj[2][3] + dMultiSilRender->viewProj[3][3]; 


		//if(pixelInfo.x<0.3f)
		{
			pixelRefine = floorf(((1.0f+fValue[0]/fValue[2])*dMultiSilRender->renderWidth)/2.0f);
			pixelRefine = min(pixelRefine, dMultiSilRender->renderWidth-1);
			pixelRefine = max(0, pixelRefine);

			//pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelRefine);
			//(*pixel++) = 0.0f;
			//(*pixel++) = 0.0f;
			//(*pixel++) = 0.0f;
			//int aaa = (patchIndexX+patchIndexY+patchIndexZ)%3;
			//pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelRefine + aaa);
			//(*pixel) = 1.0f;

			pixel = ((float*)(surfaceOut+pixelY*pitch)+ 4*pixelX + 2);
			dir = (pixelRefine>pixelX) ? 1 : -1;
			pad = *(pixelInfo-4*dir);
			//pad = 0.0f;
			//for(pixelIter=pixelX; (pixelIter!=pixelRefine+dir); pixelIter+=dir)
			//{
			//	if((*pixel)==0.0f)
			//	{
			//		(*pixel--) = 0.5f;
			//		(*pixel)   = pad;
			//	}
			//	pixel += (4*dir+1);
			//}
			int pixelPos = pixelY*dMultiSilRender->renderWidth+pixelX;
			int patchI = patchIndexX*segShift.x+patchIndexY*segShift.y+patchIndexZ*segShift.z;
			for(pixelIter=pixelX; (pixelIter!=pixelRefine+dir); pixelIter+=dir)
			{
				//if((*pixel)!=0.5f)
				{
					//(*pixel--) = 0.5f;
					(*pixel)   = pad;
				}
				dPointOccupied[pixelPos] = (((pixelIter<<pixelBits)|pixelY)<<patchBits)|patchI;
				dPointOccupiedValid[pixelPos] = 1;
				pixelPos += dir;
				pixel += (4*dir);
			}

			//(*pixel--) = 0.0f;
			//(*pixel--) = 0.0f;
			//(*pixel--) = 1.0f;
			//pixel-=4*dir;
			//(*pixel++) = 1.0f;
			//(*pixel++) = 0.0f;
			//(*pixel++) = 0.0f;

			OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000 + (pixelRefine-pixelX);
			
		}
		//else
		//{
		//	pixelRefine = floorf(((1.0f-fValue[1]/fValue[2])*dMultiSilRender->renderHeight)/2.0f);
		//	dir = (pixelRefine>pixelY) ? 1 : -1;
		//	pad = (pixelRefine>pixelY) ? pixelInfo.y : pixelInfo.w;		
		//	//pad = 0.0f;
		//	for(pixelIter=pixelY; (pixelIter!=pixelRefine&&pixelIter>=0&&pixelIter<dMultiSilRender->renderHeight); pixelIter+=dir, pixel+=(pitch/4)*dir)
		//		(*pixel) = pad;
		//	if(pixelIter==pixelRefine)
		//		(*pixel) = pad;
		//	//OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = pixelRefine-pixelY;
		//}
	}
	//else
	//{
	//	pixel = ((float*)(surfaceOut+pixelY*pitch)+ 4*pixelX);
	//	(*pixel++) = 0.0f;
	//	(*pixel++) = 0.0f;
	//	(*pixel++) = 1.0f;
	//}
	//return;


	

	//if(errF<dMultiSilRender->errorTol || deltaPoint<dMultiSilRender->errorXTol)
	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	//{
	//	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	//	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	//	xyz[2] = (xyz[2]-patchZ)/boxSize.z;
	//}
	//else
	//{
	//	xyz[0] += fValue[0];
	//	xyz[1] += fValue[0];
	//	xyz[2] += fValue[0];
	//	patchIndexX -= step[0];
	//	patchIndexY -= step[1];
	//	patchIndexZ -= step[2];
	//}

	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	patchX = patchIndexX*boxSize.x + box0.x;
	patchY = patchIndexY*boxSize.y + box0.y;
	patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;
	

	plane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	plane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	plane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	plane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*plane[0] + patchY*plane[1] + patchZ*plane[2];
	plane[0] *= boxSize.x;
	plane[1] *= boxSize.y;
	plane[2] *= boxSize.z;
	fValue[0] = fValue[1] = fValue[2] = 0.0f;
	errF = 1000.0f;
	deltaPoint = 1000.0f;

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];
		
		step[0] = floorf(xyz[0]);
		step[1] = floorf(xyz[1]);
		step[2] = floorf(xyz[2]);

		patchIndexX += step[0];
		patchIndexY += step[1];
		patchIndexZ += step[2];

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		polarIndex = PolarIndex[patchIndex];
		if(polarIndex<0)		break;

		patchIndex *= degSize;
		
		patchX += step[0]*boxSize.x;
		patchY += step[1]*boxSize.y;
		patchZ += step[2]*boxSize.z;

		plane[3] += step[0]*plane[0] + step[1]*plane[1] + step[2]*plane[2];

		xyz[0] -= step[0];
		xyz[1] -= step[1];
		xyz[2] -= step[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &fValue[1]);
		fValue[2] = plane[0]*xyz[0] + plane[1]*xyz[1] + plane[2]*xyz[2] + plane[3];


		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		EvalBezierBlossomPolarTex(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][0]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, patchIndex, &jacobi[1][1]);
		EvalBezierBlossomPolarTex(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, patchIndex, &jacobi[1][2]);

		jacobi[1][0] = (jacobi[1][0]-fValue[1])/delta;
		jacobi[1][1] = (jacobi[1][1]-fValue[1])/delta;
		jacobi[1][2] = (jacobi[1][2]-fValue[1])/delta;

		jacobi[2][0] = plane[0];
		jacobi[2][1] = plane[1];
		jacobi[2][2] = plane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= ERROR_FUNC)					break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= ERROR_DELTA)			break;	
	}
	

	isOK = true;

	//if(pixelX==14&&pixelY==133)
	//{
	//	OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000000;
	//	return;
	//}

	if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 0;
		isOK = false;
		//return;
	}

	if(patchIndex<0)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1;
		isOK = false;
		//return;
	}
	
	if(polarIndex<0)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 2;
		isOK = false;
		//return;
	}

	if(errF>ERROR_FUNC)
	{
		OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 3;
		isOK=false;
		//return;
	}

	

	//if(errF<dMultiSilRender->errorTol*5 && deltaPoint<dMultiSilRender->errorXTol*5)
	//if(errF<dMultiSilRender->errorTol)
	//if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*100 && deltaPoint<dMultiSilRender->errorXTol) )
	if(isOK)
	{
		xyz[0] = xyz[0]*boxSize.x+patchX;
		xyz[1] = xyz[1]*boxSize.y+patchY;
		xyz[2] = xyz[2]*boxSize.z+patchZ;

		//fValue[0] = xyz[0]*dMultiSilRender->viewProj[0][0] + xyz[1]*dMultiSilRender->viewProj[1][0] + xyz[2]*dMultiSilRender->viewProj[2][0] + dMultiSilRender->viewProj[3][0]; 
		fValue[1] = xyz[0]*dMultiSilRender->viewProj[0][1] + xyz[1]*dMultiSilRender->viewProj[1][1] + xyz[2]*dMultiSilRender->viewProj[2][1] + dMultiSilRender->viewProj[3][1]; 
		fValue[2] = xyz[0]*dMultiSilRender->viewProj[0][3] + xyz[1]*dMultiSilRender->viewProj[1][3] + xyz[2]*dMultiSilRender->viewProj[2][3] + dMultiSilRender->viewProj[3][3]; 

		//if(pixelInfo.x<0.3f)
		//{
		//	pixelRefine = floorf(((1.0f+fValue[0]/fValue[2])*dMultiSilRender->renderWidth)/2.0f);
		//	dir = (pixelRefine>pixelX) ? 1 : -1;
		//	pad = (pixelRefine>pixelX) ? *(pixelInfo-2) : *(pixelInfo+6);
		//	//pad = 1.0f;
		//	for(pixelIter=pixelX; (pixelIter!=pixelRefine&&pixelIter>=0&&pixelIter<dMultiSilRender->renderWidth); pixelIter+=dir, pixel+=4*dir)
		//		(*pixel) = pad;
		//	if(pixelIter==pixelRefine)
		//		(*pixel) = pad;
		//	OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = 1000 + (pixelRefine-pixelX);
		//}
		//else
		{
			pixelRefine = floorf(((1.0f-fValue[1]/fValue[2])*dMultiSilRender->renderHeight)/2.0f);
			pixelRefine = min(pixelRefine, dMultiSilRender->renderHeight-1);
			pixelRefine = max(0, pixelRefine);
			
			pixel = ((float*)(surfaceOut+pixelY*pitch)+4*pixelX+2);
			dir = (pixelRefine>pixelY) ? 1 : -1;
			//if(pixelY>0)
			//	pad = *((float*)(surfaceIn+(pixelY-dir)*pitch)+4*pixelX+2);
			//else
				pad = *((float*)(surfaceOut+pixelY*pitch)+4*pixelX+2);
			//pad = *(pixelInfo-pitch/4*dir);		
			//pad = 0.0f;

			int pixelPos = pixelY*dMultiSilRender->renderWidth+pixelX;
			int patchI = patchIndexX*segShift.x+patchIndexY*segShift.y+patchIndexZ*segShift.z;

			for(pixelIter=pixelY; (pixelIter!=pixelRefine+dir); pixelIter+=dir, pixel+=(pitch/4)*dir)
			{
				dPointOccupied[pixelPos] = (((pixelX<<pixelBits)|pixelIter)<<patchBits)|patchI;
				dPointOccupiedValid[pixelPos] = 1;
				pixelPos += dMultiSilRender->renderWidth*dir;
				(*pixel) = pad;
			}
			//(*pixel) = *((float*)(surfaceIn+pixelY*pitch)+4*pixelX+2);

			//(*pixel--) = 0.0f;
			//(*pixel--) = 0.0f;
			//(*pixel--) = 1.0f;
			//pixel-=(pitch/4)*dir;
			//(*pixel++) = 1.0f;
			//(*pixel++) = 0.0f;
			//(*pixel++) = 0.0f;

			OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x] = pixelRefine;
		}
	}
	//else
	//{
	//	pixel = ((float*)(surfaceOut+pixelY*pitch)+ 4*pixelX);
	//	(*pixel++) = 0.0f;
	//	(*pixel++) = 0.0f;
	//	(*pixel++) = 1.0f;
	//}

/**/	
}

__global__ void NoSortShadingKernel(//unsigned char*surfaceIn, 
									uint* debugInfo,
									unsigned char*surfaceOut, size_t pitch,
									int renderWidth, int renderHeight,
									MULTI_SIL_RENDER* dMultiSilRender, uint pixelBits,
									float* FuncData, int* PatchIndex,
									uint3 deg, uint degSize, uint3 seg, uint3 segShift,
									float3 boxSize, float3 box0)
{

	int pixelX = blockIdx.x*blockDim.x + threadIdx.x;
	int pixelY = blockIdx.y*blockDim.y + threadIdx.y;

	if(pixelX>=renderWidth || pixelY>=renderHeight)		return;

	float* pixelInfo = ((float*)(surfaceOut+pixelY*pitch) + pixelX*4);
	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *(pixelInfo+2);

	if(screenZ <= 0.0001f )									
	{
		pixelInfo[0] = 0;
		pixelInfo[1] = 0;
		pixelInfo[2] = 0;
		return;
	}


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float horPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float verPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

	float worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;
	
	int patchIndexX = floorf((xyz[0]-box0.x)/boxSize.x);
	int patchIndexY = floorf((xyz[1]-box0.y)/boxSize.y);
	int patchIndexZ = floorf((xyz[2]-box0.z)/boxSize.z);

	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;

	horPlane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
	horPlane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
	horPlane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
	horPlane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*horPlane[0] + patchY*horPlane[1] + patchZ*horPlane[2];
	horPlane[0] *= boxSize.x;
	horPlane[1] *= boxSize.y;
	horPlane[2] *= boxSize.z;
	verPlane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	verPlane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	verPlane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	verPlane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*verPlane[0] + patchY*verPlane[1] + patchZ*verPlane[2];
	verPlane[0] *= boxSize.x;
	verPlane[1] *= boxSize.y;
	verPlane[2] *= boxSize.z;

	int patchIndex;
	float errF = 10.0f;
	float deltaPoint = 0.0f;
	float delta = 0.01f;
	int step[3] = {0, 0, 0};

	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		step[0] = floorf(xyz[0]);
		step[1] = floorf(xyz[1]);
		step[2] = floorf(xyz[2]);

		patchIndexX += step[0];
		patchIndexY += step[1];
		patchIndexZ += step[2];

		if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
			break;
		patchIndex = patchIndexX*segShift.x + patchIndexY*segShift.y + patchIndexZ*segShift.z;
		patchIndex = PatchIndex[patchIndex];
		if(patchIndex<0)		break;
		patchIndex *= degSize;
		
		patchX += step[0]*boxSize.x;
		patchY += step[1]*boxSize.y;
		patchZ += step[2]*boxSize.z;

		horPlane[3] += step[0]*horPlane[0] + step[1]*horPlane[1] + step[2]*horPlane[2];
		verPlane[3] += step[0]*verPlane[0] + step[1]*verPlane[1] + step[2]*verPlane[2];

		xyz[0] -= step[0];
		xyz[1] -= step[1];
		xyz[2] -= step[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		fValue[1] = horPlane[0]*xyz[0] + horPlane[1]*xyz[1] + horPlane[2]*xyz[2] + horPlane[3];
		fValue[2] = verPlane[0]*xyz[0] + verPlane[1]*xyz[1] + verPlane[2]*xyz[2] + verPlane[3];


		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];

		jacobi[1][0] = horPlane[0];
		jacobi[1][1] = horPlane[1];
		jacobi[1][2] = horPlane[2];
		jacobi[2][0] = verPlane[0];
		jacobi[2][1] = verPlane[1];
		jacobi[2][2] = verPlane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= ERROR_FUNC)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabsf(fValue[0])+fabsf(fValue[1])+fabsf(fValue[2]);
		if (deltaPoint <= ERROR_DELTA)			break;	
	}

	//if(pixelX==138&&pixelY==31)
	//{
	//	debugInfo[pixelY*renderWidth+pixelX] = 10000;
	//	return;
	//}
	
		//debugInfo[pixelY*renderWidth+pixelX] = int(errF*10000);
		//return;



	bool isOK = true;
	if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
	{
		debugInfo[pixelY*renderWidth+pixelX] = 0;
		isOK = false;
		//return;
	}

	if(patchIndex<0)	
	{
		debugInfo[pixelY*renderWidth+pixelX] = 1;
		isOK = false;
		//return;
	}

	if(errF>ERROR_FUNC)
	{
		debugInfo[pixelY*renderWidth+pixelX] = 2;
		isOK = false;
		//return;
	}


	//if(errF<=dMultiSilRender->errorTol*10)
	//if(errF<dMultiSilRender->errorTol || (errF<dMultiSilRender->errorTol*100 && deltaPoint<dMultiSilRender->errorXTol) )
	//if(isOK&&errF<dMultiSilRender->errorTol*1000000)
	if(isOK)
	{
	MATERIAL material;
	material.ka = dMultiSilRender->material.ka;
	material.kd = dMultiSilRender->material.kd;
	material.ks = dMultiSilRender->material.ks;
	material.a[0] = dMultiSilRender->material.a[0];
	material.a[1] = dMultiSilRender->material.a[1];
	material.a[2] = dMultiSilRender->material.a[2];
	material.s[0] = dMultiSilRender->material.s[0];
	material.s[1] = dMultiSilRender->material.s[1];
	material.s[2] = dMultiSilRender->material.s[2];

	switch ((patchIndexX+patchIndexY+patchIndexZ)%4) 
	{
		case  0:
			material.d[0]=1.0; material.d[1]=material.d[2]=0.0;
			break;
		case  1:
			material.d[1]=1.0; material.d[0]=material.d[2]=0.0;
			break;
		case  3:
			material.d[2]=1.0; material.d[1]=material.d[0]=0.0;
			break;
		case  2:
			material.d[1]=0.0; material.d[2]=material.d[0]=1.0;
			break;
	}

	float rgb[3]={0.0f, 0.0f, 0.0f};

	norm[0] /= boxSize.x;
	norm[1] /= boxSize.y;
	norm[2] /= boxSize.z;
	normalize(norm);
	xyz[0] = xyz[0]*boxSize.x+patchX;
	xyz[1] = xyz[1]*boxSize.y+patchY;
	xyz[2] = xyz[2]*boxSize.z+patchZ;
    

	float ray[3]={xyz[0]-dMultiSilRender->eye[0], xyz[1]-dMultiSilRender->eye[1], xyz[2]-dMultiSilRender->eye[2]};
	normalize(ray);

    float diff;
	for(int i=0; i<dMultiSilRender->nLight; i++)  
	{
        diff = -(norm[0]*dMultiSilRender->light[i].dir[0]+norm[1]*dMultiSilRender->light[i].dir[1]+norm[2]*dMultiSilRender->light[i].dir[2]);
        if (diff>0)
        {
            //addPhong(&dMultiSilRender->material, &dMultiSilRender->light[i], norm, ray, rgb);
            //addDiffuse(&dMultiSilRender->material, &dMultiSilRender->light[i], diff, norm, rgb);
			addDiffuse(&material, &dMultiSilRender->light[i], diff, norm, rgb);
        }
	}
	//addambient(&dMultiSilRender->material, rgb);
	addambient(&material, rgb);

	//pixelInfo = ((float*)(surfaceOut+pixelY*pitch) + pixelX*4);
    pixelInfo[0] = rgb[0];
    pixelInfo[1] = rgb[1];
    pixelInfo[2] = rgb[2];
	}
	else
	{
		pixelInfo[0] = 0;
		pixelInfo[1] = 0;
		pixelInfo[2] = 0;

	}

}

__global__ void UpdatePointKernel(
								  //unsigned char* surfaceIn, 
								  unsigned char* surfaceOut, size_t pitch, 
								  //int width, int height, float3 eye, float errorXTol, 
								  //float4 viewProj0, float4 viewProj1, float4 viewProj2, float4 viewProj3,
								  //float4 invViewProj0, float4 invViewProj1, float4 invViewProj2, float4 invViewProj3,
								  MULTI_SIL_RENDER* dMultiSilRender, //const uint allReadRenderNum, const uint otherReadRenderNum,
								  uint pixelBits, uint patchBits, 
								  //uint* SilOccupiedSort, uint*SilOccupiedSortIndex,
								  uint* PointOccupiedSortPad, uint* PointOccupied, uint* PointOccupiedValid,
								  float* FuncData, int* PatchIndex,
								  uint3 deg, uint3 degShift, uint3 seg, uint3 segShift, 
								  float3 boxSize, float3 box0, int iter)
{	
	//uint3 degShift;
	//degShift.z = 1;
	//degShift.y = deg.z+1;
	//degShift.x = (deg.y+1)*degShift.y;
	uint ind = threadIdx.x;
	__shared__ uint patchInfo[PIXELTHREADS];
	patchInfo[ind] = PointOccupiedSortPad[ind+blockIdx.x*blockDim.x];
	__syncthreads();

	uint patchIndex = (patchInfo[0]<<(pixelBits*2))>>(pixelBits*2);
	uint patchIndexBak = patchIndex;
	uint patchIndexX = patchIndex/segShift.x;
	uint patchIndexY = patchIndex%segShift.x/segShift.y;
	uint patchIndexZ = patchIndex%segShift.y;
	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	patchIndex = PatchIndex[patchIndex];

	extern __shared__ float funcNormal[];
	uint funcSize = (deg.x+1)*degShift.x;
	float* func = funcNormal;
#ifdef NORMALFUNCTION
	float* funcJacobi = &func[funcSize];
#else
#endif

	if(ind<funcSize)
		func[ind] = FuncData[patchIndex*funcSize+ind];
	__syncthreads();

#ifdef NORMALFUNCTION
	if((ind/degShift.x)<deg.x)
	{
		funcJacobi[ind]				= deg.x*(func[ind +degShift.x] - func[ind] );
	}
	if((ind%degShift.x/degShift.y)<deg.y)
	{
		funcJacobi[funcSize +ind]	= deg.y*(func[ind +degShift.y] - func[ind] );
	}
	if((ind%degShift.y)<deg.z)
	{
		funcJacobi[funcSize*2 +ind] = deg.z*(func[ind +degShift.z] - func[ind] );
	}
	__syncthreads();
#else
	float delta = 0.01f;
#endif

	if(patchInfo[ind]==0)
		return;


	float xyz[3] = {0.0f, 0.0f, 0.0f};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float horPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float verPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};

	uint pixelX = (patchInfo[ind]>>(patchBits+pixelBits));
	uint pixelY = (patchInfo[ind]<<pixelBits)>>(patchBits+pixelBits);

	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float screenZ = *((float*)(surfaceOut+pixelY*pitch) + pixelX*4 + 2);

	horPlane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
	horPlane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
	horPlane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
	horPlane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*horPlane[0] + patchY*horPlane[1] + patchZ*horPlane[2];
	horPlane[0] *= boxSize.x;
	horPlane[1] *= boxSize.y;
	horPlane[2] *= boxSize.z;
	verPlane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	verPlane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	verPlane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	verPlane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*verPlane[0] + patchY*verPlane[1] + patchZ*verPlane[2];
	verPlane[0] *= boxSize.x;
	verPlane[1] *= boxSize.y;
	verPlane[2] *= boxSize.z;

	float worldW = screenX*dMultiSilRender->invViewProj[0][3] + screenY*dMultiSilRender->invViewProj[1][3] + screenZ*dMultiSilRender->invViewProj[2][3] + dMultiSilRender->invViewProj[3][3];
	xyz[0] = (screenX*dMultiSilRender->invViewProj[0][0] + screenY*dMultiSilRender->invViewProj[1][0] + screenZ*dMultiSilRender->invViewProj[2][0] + dMultiSilRender->invViewProj[3][0]) / worldW;
	xyz[1] = (screenX*dMultiSilRender->invViewProj[0][1] + screenY*dMultiSilRender->invViewProj[1][1] + screenZ*dMultiSilRender->invViewProj[2][1] + dMultiSilRender->invViewProj[3][1]) / worldW;
	xyz[2] = (screenX*dMultiSilRender->invViewProj[0][2] + screenY*dMultiSilRender->invViewProj[1][2] + screenZ*dMultiSilRender->invViewProj[2][2] + dMultiSilRender->invViewProj[3][2]) / worldW;

	xyz[0] = (xyz[0]-patchX)/boxSize.x;
	xyz[1] = (xyz[1]-patchY)/boxSize.y;
	xyz[2] = (xyz[2]-patchZ)/boxSize.z;

	float* pixel = (float*)(surfaceOut + pixelY*pitch) + 4*pixelX;

	float errF = 0.0f;
	float deltaPoint = 0.0f;
	for(int its=0; its<5; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, func, &fValue[0]);
		fValue[1] = horPlane[0]*xyz[0] + horPlane[1]*xyz[1] + horPlane[2]*xyz[2] + horPlane[3];
		fValue[2] = verPlane[0]*xyz[0] + verPlane[1]*xyz[1] + verPlane[2]*xyz[2] + verPlane[3];

#ifdef NORMALFUNCTION
		EvalBezierBlossomNormal(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, degShift, funcJacobi, &jacobi[0][0]);
#else
		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, func, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, func, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, func, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;
#endif

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];
		jacobi[1][0] = horPlane[0];
		jacobi[1][1] = horPlane[1];
		jacobi[1][2] = horPlane[2];
		jacobi[2][0] = verPlane[0];
		jacobi[2][1] = verPlane[1];
		jacobi[2][2] = verPlane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= ERROR_FUNC)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabsf(fValue[0])+fabsf(fValue[1])+fabsf(fValue[2]);
		if (deltaPoint <= ERROR_DELTA)			break;	
	}

	if(xyz[0]<0||xyz[1]<0||xyz[2]<0||xyz[0]>1||xyz[1]>1||xyz[2]>1)
	{
		int step[3];
		step[0] = floorf(xyz[0]);
		step[1] = floorf(xyz[1]);
		step[2] = floorf(xyz[2]);

		pixel[0] = xyz[0] - step[0];
		pixel[1] = xyz[1] - step[1];
		pixel[2] = xyz[2] - step[2];

		//if(patchIndexX<0 || patchIndexX>=seg.x || patchIndexY<0 || patchIndexY>=seg.y || patchIndexZ<0 || patchIndexZ>=seg.z)
		int patchI = patchIndexBak + step[0]*segShift.x + step[1]*segShift.y + step[2]*segShift.z;
		int pixelI = pixelY*dMultiSilRender->renderWidth+pixelX;
		//if(patchI>0)
		{
			PointOccupied[pixelI] = (((pixelX<<pixelBits)|pixelY)<<patchBits)|patchI;
			PointOccupiedValid[pixelI] = 1;
		}
		return;
	}

	MATERIAL material;
	material.ka = dMultiSilRender->material.ka;
	material.kd = dMultiSilRender->material.kd;
	material.ks = dMultiSilRender->material.ks;
	material.a[0] = dMultiSilRender->material.a[0];
	material.a[1] = dMultiSilRender->material.a[1];
	material.a[2] = dMultiSilRender->material.a[2];
	material.s[0] = dMultiSilRender->material.s[0];
	material.s[1] = dMultiSilRender->material.s[1];
	material.s[2] = dMultiSilRender->material.s[2];

	switch ((patchIndexX+patchIndexY+patchIndexZ)%4) 
	{
		case  0:
			material.d[0]=1.0; material.d[1]=material.d[2]=0.0;
			break;
		case  1:
			material.d[1]=1.0; material.d[0]=material.d[2]=0.0;
			break;
		case  3:
			material.d[2]=1.0; material.d[1]=material.d[0]=0.0;
			break;
		case  2:
			material.d[1]=0.0; material.d[2]=material.d[0]=1.0;
			break;
	}

    float rgb[3]={0.0f, 0.0f, 0.0f};

	norm[0] /= boxSize.x;
	norm[1] /= boxSize.y;
	norm[2] /= boxSize.z;
	normalize(norm);
	xyz[0] = xyz[0]*boxSize.x+patchX;
	xyz[1] = xyz[1]*boxSize.y+patchY;
	xyz[2] = xyz[2]*boxSize.z+patchZ;
    

	float ray[3]={xyz[0]-dMultiSilRender->eye[0], xyz[1]-dMultiSilRender->eye[1], xyz[2]-dMultiSilRender->eye[2]};
	normalize(ray);

    float diff;
	for(int i=0; i<dMultiSilRender->nLight; i++)  
	{
        diff = -(norm[0]*dMultiSilRender->light[i].dir[0]+norm[1]*dMultiSilRender->light[i].dir[1]+norm[2]*dMultiSilRender->light[i].dir[2]);
        if (diff>0)
        {
            //addPhong(&dMultiSilRender->material, &dMultiSilRender->light[i], norm, ray, rgb);
            //addDiffuse(&dMultiSilRender->material, &dMultiSilRender->light[i], diff, norm, rgb);
			addDiffuse(&material, &dMultiSilRender->light[i], diff, norm, rgb);
        }
	}
	//addambient(&dMultiSilRender->material, rgb);
	addambient(&material, rgb);
    pixel[0] = rgb[0];
    pixel[1] = rgb[1];
    pixel[2] = rgb[2];
}

__global__ void UpdateBoundaryPointKernel(
								  unsigned char* surfaceOut, size_t pitch, 
								  //int width, int height, float3 eye, float errorXTol, 
								  //float4 viewProj0, float4 viewProj1, float4 viewProj2, float4 viewProj3,
								  //float4 invViewProj0, float4 invViewProj1, float4 invViewProj2, float4 invViewProj3,
								  MULTI_SIL_RENDER* dMultiSilRender, 
								  uint pixelBits, uint patchBits, 
								  uint* OccupiedCompact,
								  float* FuncData, int* PatchIndex, uint nPatch,
								  uint3 deg, uint3 degShift, uint3 seg, uint3 segShift, 
								  float3 boxSize, float3 box0, int iter)
{	
	uint patchInfo = OccupiedCompact[threadIdx.x+blockIdx.x*blockDim.x];

	if(patchInfo==0)
		return;

	uint patchIndex = (patchInfo<<(pixelBits*2))>>(pixelBits*2);
	uint patchIndexX = patchIndex/segShift.x;
	uint patchIndexY = patchIndex%segShift.x/segShift.y;
	uint patchIndexZ = patchIndex%segShift.y;
	float patchX = patchIndexX*boxSize.x + box0.x;
	float patchY = patchIndexY*boxSize.y + box0.y;
	float patchZ = patchIndexZ*boxSize.z + box0.z;

	patchIndex = PatchIndex[patchIndex];
	if(patchIndex<0 || patchIndex>=nPatch)
		return;
	patchIndex *= (deg.x+1)*degShift.x;

	uint pixelX = (patchInfo>>(patchBits+pixelBits));
	uint pixelY = (patchInfo<<pixelBits)>>(patchBits+pixelBits);
	float screenX = 2.0f*(pixelX+0.5f)/dMultiSilRender->renderWidth-1.0f;
	float screenY = 1.0f-2.0f*(pixelY+0.5f)/dMultiSilRender->renderHeight;
	float* pixel = (float*)(surfaceOut + pixelY*pitch) + 4*pixelX;
	float xyz[3] = {pixel[0], pixel[1], pixel[2]};
	float norm[3] = {0.0f, 0.0f, 0.0f};
	float fValue[3] = {0.0f, 0.0f, 0.0f};
	float horPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float verPlane[4] = {0.0f, 0.0f, 0.0f, 0.0f};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	float errF = 1000.0f;
	float deltaPoint = 1000.0f;
	float delta = 0.01f;

	horPlane[0] = (screenY*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][1]);
	horPlane[1] = (screenY*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][1]);
	horPlane[2] = (screenY*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][1]);
	horPlane[3] = (screenY*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][1]) + patchX*horPlane[0] + patchY*horPlane[1] + patchZ*horPlane[2];
	horPlane[0] *= boxSize.x;
	horPlane[1] *= boxSize.y;
	horPlane[2] *= boxSize.z;
	verPlane[0] = (screenX*dMultiSilRender->viewProj[0][3]-dMultiSilRender->viewProj[0][0]);
	verPlane[1] = (screenX*dMultiSilRender->viewProj[1][3]-dMultiSilRender->viewProj[1][0]);
	verPlane[2] = (screenX*dMultiSilRender->viewProj[2][3]-dMultiSilRender->viewProj[2][0]);
	verPlane[3] = (screenX*dMultiSilRender->viewProj[3][3]-dMultiSilRender->viewProj[3][0]) + patchX*verPlane[0] + patchY*verPlane[1] + patchZ*verPlane[2];
	verPlane[0] *= boxSize.x;
	verPlane[1] *= boxSize.y;
	verPlane[2] *= boxSize.z;

	for(int its=0; its<3; its++)
	{
		xyz[0] -= fValue[0];
		xyz[1] -= fValue[1];
		xyz[2] -= fValue[2];

		EvalBezierBlossom(xyz[0], xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &fValue[0]);
		fValue[1] = horPlane[0]*xyz[0] + horPlane[1]*xyz[1] + horPlane[2]*xyz[2] + horPlane[3];
		fValue[2] = verPlane[0]*xyz[0] + verPlane[1]*xyz[1] + verPlane[2]*xyz[2] + verPlane[3];

		EvalBezierBlossom(xyz[0]+delta, xyz[1], xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][0]);
		EvalBezierBlossom(xyz[0], xyz[1]+delta, xyz[2], deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][1]);
		EvalBezierBlossom(xyz[0], xyz[1], xyz[2]+delta, deg.x, deg.y, deg.z, FuncData+patchIndex, &jacobi[0][2]);
		jacobi[0][0] = (jacobi[0][0]-fValue[0])/delta;
		jacobi[0][1] = (jacobi[0][1]-fValue[0])/delta;
		jacobi[0][2] = (jacobi[0][2]-fValue[0])/delta;

		norm[0] = jacobi[0][0];
		norm[1] = jacobi[0][1];
		norm[2] = jacobi[0][2];

		jacobi[1][0] = horPlane[0];
		jacobi[1][1] = horPlane[1];
		jacobi[1][2] = horPlane[2];
		jacobi[2][0] = verPlane[0];
		jacobi[2][1] = verPlane[1];
		jacobi[2][2] = verPlane[2];

		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= ERROR_FUNC)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabsf(fValue[0])+fabsf(fValue[1])+fabsf(fValue[2]);
		if (deltaPoint <= ERROR_DELTA)			break;	
	}

	

	MATERIAL material;
	material.ka = dMultiSilRender->material.ka;
	material.kd = dMultiSilRender->material.kd;
	material.ks = dMultiSilRender->material.ks;
	material.a[0] = dMultiSilRender->material.a[0];
	material.a[1] = dMultiSilRender->material.a[1];
	material.a[2] = dMultiSilRender->material.a[2];
	material.s[0] = dMultiSilRender->material.s[0];
	material.s[1] = dMultiSilRender->material.s[1];
	material.s[2] = dMultiSilRender->material.s[2];

	switch ((patchIndexX+patchIndexY+patchIndexZ)%4) 
	{
		case  0:
			material.d[0]=1.0; material.d[1]=material.d[2]=0.0;
			break;
		case  1:
			material.d[1]=1.0; material.d[0]=material.d[2]=0.0;
			break;
		case  3:
			material.d[2]=1.0; material.d[1]=material.d[0]=0.0;
			break;
		case  2:
			material.d[1]=0.0; material.d[2]=material.d[0]=1.0;
			break;
	}

    float rgb[3]={0.0f, 0.0f, 0.0f};

	norm[0] /= boxSize.x;
	norm[1] /= boxSize.y;
	norm[2] /= boxSize.z;
	normalize(norm);
	xyz[0] = xyz[0]*boxSize.x+patchX;
	xyz[1] = xyz[1]*boxSize.y+patchY;
	xyz[2] = xyz[2]*boxSize.z+patchZ;
    

	float ray[3]={xyz[0]-dMultiSilRender->eye[0], xyz[1]-dMultiSilRender->eye[1], xyz[2]-dMultiSilRender->eye[2]};
	normalize(ray);

    float diff;
	for(int i=0; i<dMultiSilRender->nLight; i++)  
	{
        diff = -(norm[0]*dMultiSilRender->light[i].dir[0]+norm[1]*dMultiSilRender->light[i].dir[1]+norm[2]*dMultiSilRender->light[i].dir[2]);
        if (diff>0)
        {
            //addPhong(&dMultiSilRender->material, &dMultiSilRender->light[i], norm, ray, rgb);
            //addDiffuse(&dMultiSilRender->material, &dMultiSilRender->light[i], diff, norm, rgb);
			addDiffuse(&material, &dMultiSilRender->light[i], diff, norm, rgb);
        }
	}
	//addambient(&dMultiSilRender->material, rgb);
	addambient(&material, rgb);
    pixel[0] = rgb[0];
    pixel[1] = rgb[1];
    pixel[2] = rgb[2];
}




extern "C"
void MultiSilRender(MULTI_SIL_RENDER* multiSilRender, void* surfaceIn, void* surfaceOut, size_t pitch, int iter)
{
	CUDA_SAFE_CALL( cudaMemcpy(dMultiSilRender, multiSilRender, sizeof(MULTI_SIL_RENDER), cudaMemcpyHostToDevice) );
	invViewProj0.x = multiSilRender->invViewProj[0][0]; invViewProj0.y = multiSilRender->invViewProj[0][1]; invViewProj0.z = multiSilRender->invViewProj[0][2]; invViewProj0.w = multiSilRender->invViewProj[0][3]; 
	invViewProj1.x = multiSilRender->invViewProj[1][0]; invViewProj1.y = multiSilRender->invViewProj[1][1]; invViewProj1.z = multiSilRender->invViewProj[1][2]; invViewProj1.w = multiSilRender->invViewProj[1][3]; 
	invViewProj2.x = multiSilRender->invViewProj[2][0]; invViewProj2.y = multiSilRender->invViewProj[2][1]; invViewProj2.z = multiSilRender->invViewProj[2][2]; invViewProj2.w = multiSilRender->invViewProj[2][3]; 
	invViewProj3.x = multiSilRender->invViewProj[3][0]; invViewProj3.y = multiSilRender->invViewProj[3][1]; invViewProj3.z = multiSilRender->invViewProj[3][2]; invViewProj3.w = multiSilRender->invViewProj[3][3]; 


	CUDA_SAFE_CALL( cudaMemset(dPolarIndex, -1, nPatch*sizeof(int)) );

	float3 eye = make_float3(multiSilRender->eye[0], multiSilRender->eye[1], multiSilRender->eye[2]);
	size_t Ns = degSize*20*sizeof(float);
#if	UNIFORM_KNOT
	//GetPolarData<<<segSize, degSize, Ns>>>(dFuncData, dPolarData, dPatchIndex, dPolarIndex, eye, deg, degShift, segShift, degSize, box0, boxSize);
	GetPolarData<<<nPatch, degSize, Ns>>>(dFuncData, dPolarData, /*dPatchIndex, */dPatchIndexInverse, dPolarIndex, eye, deg, degShift, segShift, degSize, box0, boxSize);
#else
#endif


#if DEBUG_BUFFERS
	printf("dPolarData:\n");
	dumpBufferfloat(dPolarData, degSize*nPatch);
#endif

	//int silSize = renderWidth*renderHeight*sizeof(uint);
	CUDA_SAFE_CALL( cudaMemset(dSilOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dSilOccupiedValid, 0, allPointSize) );
	dim3 Db = dim3(16,16);
	dim3 Dg = dim3( renderWidth/Db.x+1, renderHeight/Db.y+1 );

#if	UNIFORM_KNOT	
#if SORT_SHARE
	CUDA_SAFE_CALL( cudaMemset(dPointOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dPointOccupiedValid, 0, allPointSize) );
	CheckAllSilPointUniform<<<Dg,Db>>>((unsigned char*)surfaceIn, pitch, renderWidth, renderHeight, 
								pixelBits, dSilOccupied, dSilOccupiedValid,
								box0, boxSize, seg, segShift, patchBits,
								dPatchIndex, dPointOccupied, dPointOccupiedValid, 
								invViewProj0, invViewProj1, invViewProj2, invViewProj3);
#else
	NoSortCheckSilPointUniform<<<Dg,Db>>>((unsigned char*)surfaceIn, pitch, renderWidth, renderHeight, 
									pixelBits, dSilOccupied, dSilOccupiedValid);	
#endif
#else
#endif

#if DEBUG_BUFFERS
	printf("dSilOccupied:\n");
	dumpBuffer(dPointOccupiedValid, renderWidth*renderHeight);
#endif

	CUDA_SAFE_CALL( cudaMemset(dOccupiedCompact, 0, allPointSize) );

	CUDPPConfiguration config;
	CUDPPResult result = CUDPP_SUCCESS;	

    config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, renderWidth*renderHeight, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, dOccupiedCompact, dMaxSilPoint, dSilOccupied, dSilOccupiedValid, renderWidth*renderHeight);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }

	actualSilPoint=0;
	CUDA_SAFE_CALL( cudaMemcpy(&actualSilPoint, dMaxSilPoint, sizeof(size_t), cudaMemcpyDeviceToHost) );

#if DEBUG_BUFFERS
	float* elimation;
	CUDA_SAFE_CALL( cudaMalloc((void**)&elimation, actualSilPoint*15*sizeof(float)) );
	printf("dOccupiedCompact:\n");
	dumpBuffer(dOccupiedCompact, actualSilPoint);
#endif

#if	UNIFORM_KNOT	
	NoSortUpdateSilPointKernel<<<actualSilPoint/PIXELTHREADS+1, PIXELTHREADS>>>
								(dPointOccupied, dPointOccupiedValid, patchBits,
								(unsigned char*)surfaceIn, (unsigned char*)surfaceOut, pitch,
								dMultiSilRender, pixelBits,
								dOccupiedCompact, dFuncData, dPatchIndex, dPolarIndex,
								deg, degSize, seg, segShift,
								boxSize, box0//,elimation
								);
#else
#endif

#if DEBUG_BUFFERS
	printf("dOccupiedCompact:\n");
	dumpBufferfloat(elimation, actualSilPoint*15);
	CUDA_SAFE_CALL( cudaFree(elimation) );
	dumpBuffer(dOccupiedCompact, actualSilPoint);
#endif




#if	UNIFORM_KNOT	
#if SORT_SHARE

	CUDA_SAFE_CALL( cudaMemset(dOccupiedCompact, 0, allPointSize) );

    config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, renderWidth*renderHeight, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, dOccupiedCompact, dMaxPoint, dPointOccupied, dPointOccupiedValid, renderWidth*renderHeight);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}	
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }
	
	
	CUDA_SAFE_CALL( cudaMemcpy(&maxPoint, dMaxPoint, sizeof(size_t), cudaMemcpyDeviceToHost) );
#if DEBUG_BUFFERS
	printf("dOccupiedCompact:\n");
	dumpBuffer(dOccupiedCompact, maxPoint);
#endif



	CUDA_SAFE_CALL( cudaMemset(dOccupiedSort, 0, maxPoint*sizeof(uint)) );

	config.algorithm = CUDPP_SORT_RADIX_GLOBAL;
	config.datatype = CUDPP_UINT;
	config.op = CUDPP_ADD;
	config.options = (CUDPPOption)0;
	result = cudppPlan(&plan, config, maxPoint, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for SortBit\n");
		return;
	}
	result = cudppSortBit(plan, dOccupiedSort, dOccupiedCompact, maxPoint, patchBits);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error SortBit\n");
		return;
	}
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
    }


#if DEBUG_BUFFERS
	printf("dOccupiedSort:\n");
	dumpBuffer(dOccupiedSort, maxPoint);
#endif

	CUDA_SAFE_CALL( cudaMemcpy(pointOccupiedSort, dOccupiedSort, maxPoint*sizeof(uint), cudaMemcpyDeviceToHost) );
	pointOccupiedSort[maxPoint] = 0;
	uint* p1=pointOccupiedSort, *p2=pointOccupiedPad;
	int polarPad=0, local=0, np = (*p1)&patchMask;
	while(*p2++ = *p1++)
	{
		local++;
		if( (*p1&patchMask)!=np || local==PIXELTHREADS)
		{
			for(int i=local; i<PIXELTHREADS; i++)
				*p2++ = 0;
			polarPad++;
			local = 0;
			np = (*p1&patchMask);
		}		
	}
	if(local)
	{
		for(int i=local; i<PIXELTHREADS; i++)
			*p2++ = 0;
	}
	polarPad *= PIXELTHREADS;
	//printf("%d %d\n", silSize, polarPad);

	if(polarPad*sizeof(uint)>allPointSize*5/3)
		printf("Allocate more GPU memeory for dSilOccupiedSortPad!\n");

 
	CUDA_SAFE_CALL( cudaMemset(dOccupiedSortPad, 0, polarPad*sizeof(uint)) );
	CUDA_SAFE_CALL( cudaMemcpy(dOccupiedSortPad, pointOccupiedPad, polarPad*sizeof(uint), cudaMemcpyHostToDevice) );

#if DEBUG_BUFFERS
	printf("dPointOccupiedSortPad:\n");
	printf("polarPad=%d\n", polarPad);
	dumpBuffer(dOccupiedSortPad, polarPad);
#endif

	Db = dim3(PIXELTHREADS, 1, 1);
	Dg = dim3(polarPad/PIXELTHREADS, 1, 1);
	//Dg = dim3(iter+1, 1, 1);
	//Dg = dim3(nPatchEveryBlock*maxPolarIndexBlock, 1, 1);

#ifdef NORMALFUNCTION
	Ns = PIXELTHREADS*sizeof(float)*2;
#else
	//512&1024 work fine with deg.xyz=3;
	Ns =1024; //(deg.x+1)*degShift.x*3*sizeof(float)*2
#endif

	CUDA_SAFE_CALL( cudaMemset(dPointOccupied, 0, allPointSize) );
	CUDA_SAFE_CALL( cudaMemset(dPointOccupiedValid, 0, allPointSize) );

	UpdatePointKernel<<<Dg, Db, Ns>>>(
		//(unsigned char*)surfaceIn, 
		(unsigned char*)surfaceOut, pitch,
		//renderWidth, renderHeight, eye, errorXTol,
		//viewProj0, viewProj1, viewProj2, viewProj3,
		//invViewProj0, invViewProj1, invViewProj2, invViewProj3,
		dMultiSilRender, //sizeof(MULTI_SIL_RENDER)/(PIXELTHREADS*sizeof(float)), sizeof(MULTI_SIL_RENDER)%(PIXELTHREADS*sizeof(float)),
		pixelBits, patchBits, 
		//dSilOccupiedSort, dSilOccupiedSortIndex,
		dOccupiedSortPad, dPointOccupied, dPointOccupiedValid,
		dFuncData, dPatchIndex,
		deg, degShift, seg, segShift,
		boxSize, box0, iter);


	CUDA_SAFE_CALL( cudaMemset(dOccupiedCompact, 0, allPointSize) );

    config.algorithm = CUDPP_COMPACT;
    config.datatype = CUDPP_UINT;
	config.options = CUDPP_OPTION_FORWARD;
	result = cudppPlan(&plan, config, renderWidth*renderHeight, 1, 0);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error creating plan for compact\n");
		return;
	}
	result = cudppCompact(plan, dOccupiedCompact, dMaxPoint, dPointOccupied, dPointOccupiedValid, renderWidth*renderHeight);
	if (result != CUDPP_SUCCESS)
	{
		printf("Error compacting !\n");
		return;
	}	
	result = cudppDestroyPlan(plan);    
    if (result != CUDPP_SUCCESS)
    {
        printf("Error destroying CUDPPPlan!\n");
		return;
    }
	
	
	CUDA_SAFE_CALL( cudaMemcpy(&maxPoint, dMaxPoint, sizeof(size_t), cudaMemcpyDeviceToHost) );
#if DEBUG_BUFFERS
	printf("dOccupiedCompact:\n");
	dumpBuffer(dOccupiedCompact, maxPoint);
#endif

	UpdateBoundaryPointKernel<<<maxPoint/PIXELTHREADS+1, PIXELTHREADS>>>(
		(unsigned char*)surfaceOut, pitch,
		//renderWidth, renderHeight, eye, errorXTol,
		//viewProj0, viewProj1, viewProj2, viewProj3,
		//invViewProj0, invViewProj1, invViewProj2, invViewProj3,
		dMultiSilRender,
		pixelBits, patchBits, 
		dOccupiedCompact,
		dFuncData, dPatchIndex, nPatch,
		deg, degShift, seg, segShift,
		boxSize, box0, iter);

#if DEBUG_BUFFERS
	printf("dPointOccupiedValid:\n");
	dumpBuffer(dOccupiedCompact, maxPoint);
#endif


#else
	//dim3 Db, Dg;
	Db.x = Db.y = 8;
	Dg.x = renderWidth/Db.x+1;
	Dg.y = renderHeight/Db.y+1;
	NoSortShadingKernel<<<Dg, Db>>>(//(unsigned char*)surfaceIn,
									dSilOccupied,
									(unsigned char*)surfaceOut, pitch,
									renderWidth, renderHeight,
									dMultiSilRender, pixelBits,
									dFuncData, dPatchIndex,
									deg, degSize, seg, segShift,
									boxSize, box0);
#if DEBUG_BUFFERS
	printf("Shading:\n");
	dumpBuffer(dSilOccupied, renderWidth*renderHeight);
#endif

#endif

#else	//else of uinform
#endif

/**/


}

extern "C"
void CleanUpMultiRenderCUDA()
{
	CUDA_SAFE_CALL( cudaFree(dFuncData) );
	CUDA_SAFE_CALL( cudaFree(dPatchIndex) );
	CUDA_SAFE_CALL( cudaFree(dPatchIndexInverse) );

	CUDA_SAFE_CALL( cudaUnbindTexture(polarTex) );
	CUDA_SAFE_CALL( cudaFree(dPolarData) );
	CUDA_SAFE_CALL( cudaFree(dPolarIndex) );
	CUDA_SAFE_CALL( cudaFree(dMultiSilRender) );

	CUDA_SAFE_CALL( cudaFree(dSilOccupiedValid) );
	CUDA_SAFE_CALL( cudaFree(dSilOccupied) );

	CUDA_SAFE_CALL( cudaFree(dPointOccupied) );
	CUDA_SAFE_CALL( cudaFree(dPointOccupiedValid) );

	CUDA_SAFE_CALL( cudaFree(dOccupiedCompact) );
	CUDA_SAFE_CALL( cudaFree(dOccupiedSort) );

	CUDA_SAFE_CALL( cudaFree(dMaxSilPoint) );
	CUDA_SAFE_CALL( cudaFree(dMaxPoint) );

	CUDA_SAFE_CALL( cudaFree(dOccupiedSortPad) );

	free(pointOccupiedSort);
	free(pointOccupiedPad);

#ifndef UNIFORM_KNOT	
	CUDA_SAFE_CALL( cudaFree(dBoxX) );
	CUDA_SAFE_CALL( cudaFree(dBoxY) );
	CUDA_SAFE_CALL( cudaFree(dBoxZ) );
#endif
}
